#=============================================================================
# Manuscript: Pushing boundaries: How lawmakers shape judicial decision-making 
#=============================================================================

# Journal: Comparative Political Studies
# Author: Philipp Schroeder (LMU Munich)
# For questions contact P.Schroeder@gsi.uni-muenchen.de
# Date: 2021-10-08
# R version 4.0.4 (2021-02-15) -- "Lost Library Book"

### Main analyses ###

rm(list = ls())

library(rstanarm) # Version 2.21.1
library(scales)
library(texreg)

#================================================================
# A) Data

# Set your path to the folder with replication files
getwd()

source("data_processing.R") # sources data and descriptive statistics

#================================================================
# B) Estimation
options(mc.cores = parallel::detectCores())
names(mydata)

### Model 1
fit.1 <- stan_glmer(court_decision ~ obj_gov + obj_opp + auth_gov + cross_party + policy_area + (1 | law_id),
                    data = mydata, family = binomial(link = "logit"),
                    # priors (default November 2018)
                    # prior = normal(0,2.5), prior_intercept = normal(0,10),
                    chains=4, seed=123, iter=6000, warmup=1000)
print(fit.1)
#summary(fit.1)
# priors
prior_summary(fit.1)
# convergence
rhat.1 <- summary(fit.1)[, "Rhat"]
max(rhat.1)

# Sigma
var(unlist(ranef(fit.1)))
# Get full results for fixed effect coefficients
BETA.1 <- as.matrix(fit.1)[,2:length(fixef(fit.1))]
mean.BETA.1 <- apply(BETA.1,2,mean)
sd.BETA.1 <- apply(BETA.1,2,sd)
low.BETA.1 <- apply(BETA.1,2,quantile,c(0.025,0.975))[1,]
high.BETA.1 <- apply(BETA.1,2,quantile,c(0.025,0.975))[2,]
round(mean.BETA.1, digits = 2)
round(low.BETA.1, digits = 2)
round(high.BETA.1, digits = 2)



### Model 2
fit.2 <- stan_glmer(court_decision ~ obj_gov + obj_opp + z.days_passed + cross_party + policy_area + (1 | law_id),
                    data = mydata, family = binomial(link = "logit"),
                    # priors (default November 2018)
                    # prior = normal(0,2.5), prior_intercept = normal(0,10),
                    chains=4, seed=123, iter=6000, warmup=1000)
print(fit.2)
#summary(fit.2)
# priors
prior_summary(fit.2)
# convergence
rhat.2 <- summary(fit.2)[, "Rhat"]
max(rhat.2)

# Sigma
var(unlist(ranef(fit.2)))
# Get full results for fixed effect coefficients
BETA.2 <- as.matrix(fit.2)[,2:length(fixef(fit.2))]
mean.BETA.2 <- apply(BETA.2,2,mean)
sd.BETA.2 <- apply(BETA.2,2,sd)
low.BETA.2 <- apply(BETA.2,2,quantile,c(0.025,0.975))[1,]
high.BETA.2 <- apply(BETA.2,2,quantile,c(0.025,0.975))[2,]
round(mean.BETA.2, digits = 2)
round(low.BETA.2, digits = 2)
round(high.BETA.2, digits = 2)


#================================================================
# C) Predicted probabilities
names(mydata)

fixef(fit.2)
tmpdata <- as.data.frame(cbind(1, mydata$obj_gov, mydata$obj_opp, mydata$z.days_passed, mydata$cross_party,
                                  ifelse(mydata$policy_area == "education/research", 1,0),
                                  ifelse(mydata$policy_area == "family", 1,0), ifelse(mydata$policy_area == "healthcare", 1,0), ifelse(mydata$policy_area == "interior", 1,0),
                                  ifelse(mydata$policy_area == "labour/social insurance", 1,0), ifelse(mydata$policy_area == "other", 1,0), ifelse(mydata$policy_area == "public finances", 1,0),
                                  ifelse(mydata$policy_area == "rights", 1,0), ifelse(mydata$policy_area == "transport/public infrastructure", 1,0)))
dim(tmpdata)
colnames(tmpdata) <- names(fixef(fit.2))

# Set up design matrix Z for random effects
N <- dim(mydata)[1]
q <- dim(ranef(fit.2)[[1]])[1]
Z <- as.data.frame(matrix(NA, nrow = N, ncol = q))
dim(Z)
colnames(Z) <- rownames(ranef(fit.2)[[1]])
for (i in colnames(Z)) {
  Z[,i] <- ifelse(mydata$law_id == i,1,0) 
}

# Check whether data matrix and design matrix are set up correctly
mydata[1,c("case","law_id", "obj_gov","obj_opp","cross_party","policy_area")]
tmpdata[1,] # values for first observation should match mydata
Z[1,as.character(mydata[1,"law_id"])] # should return 1 if correct

# Fixed effect coefficients from all sampling iterations
BETA <- as.matrix(fit.2)[,1:length(fixef(fit.2))]
dim(BETA)
BETA[1,]

# Random effects from all sampling iterations
start.alpha <- length(fixef(fit.2)) + 1
GAMMA <- as.matrix(fit.2)[,start.alpha:dim(as.matrix(fit.2))[2]]
GAMMA <- GAMMA[,-dim(GAMMA)[2]] # drop sigma
dim(GAMMA)

# Set up data matrices with majMP.object = 0 and majMP.object = 1 for all observations
X.0 <- as.matrix(tmpdata)
X.0[,match("obj_gov",names(fixef(fit.2)))] <- 0
X.0[,match("obj_opp",names(fixef(fit.2)))] <- 0

X.1 <- as.matrix(tmpdata)
X.1[,match("obj_gov",names(fixef(fit.2)))] <- 0
X.1[,match("obj_opp",names(fixef(fit.2)))] <- 1

X.2 <- as.matrix(tmpdata)
X.2[,match("obj_gov",names(fixef(fit.2)))] <- 1
X.2[,match("obj_opp",names(fixef(fit.2)))] <- 1

Z <- as.matrix(Z)

linpred_0 <- X.0 %*% t(BETA) + Z %*% t(GAMMA) # Linear predictor with X.0
dim(linpred_0)
pred_0 <- exp(linpred_0)/(1+exp(linpred_0)) # Through link function
dim(pred_0)
avg_pred_0 <- apply(pred_0,2,mean)
mean(avg_pred_0)
quantile(avg_pred_0,c(0.025,0.975))

linpred_1 <- X.1 %*% t(BETA) + Z %*% t(GAMMA) # Linear predictor with X.1
dim(linpred_1)
pred_1 <- exp(linpred_1)/(1+exp(linpred_1)) # Through link function
dim(pred_1)
avg_pred_1 <- apply(pred_1,2,mean)
mean(avg_pred_1)
quantile(avg_pred_1,c(0.025,0.975))

linpred_2 <- X.2 %*% t(BETA) + Z %*% t(GAMMA) # Linear predictor with X.1
dim(linpred_2)
pred_2 <- exp(linpred_2)/(1+exp(linpred_2)) # Through link function
dim(pred_2)
avg_pred_2 <- apply(pred_2,2,mean)
mean(avg_pred_2)
quantile(avg_pred_2,c(0.025,0.975))

# Average prediction across observations per sampling iteration
delta_1 <- pred_1 - pred_0
dim(delta_1)
pred_delta_1 <- apply(delta_1,2,mean)
mean(pred_delta_1)
quantile(pred_delta_1,c(.025,.975))

delta_2 <- pred_2 - pred_1
dim(delta_2)
pred_delta_2 <- apply(delta_2,2,mean)
mean(pred_delta_2)
quantile(pred_delta_2,c(.025,.975))

pdf("probabilities.pdf", height = 5, width = 7) # Figure 5
par(mar=c(7,4,4,2))
plot("",
     xlim = c(0.75,3),
     ylim = c(0,1),
     xlab = "",
     ylab = "Predicted probability",
     main = "Predicted probabilities of Strike = 1",
     axes = FALSE)
xtick <- seq(1, 3, 1)
axis(side = 1, at = xtick, labels = FALSE)
text(x=xtick, y = -0.075, 
     labels = c("No objections","Objections by \n opposition lawmakers","Objections by opposition \n and governing lawmakers"), srt = 45, pos = 2, xpd = TRUE, cex = 0.8)
axis(2,at = seq(0,1, by = .2))
segments(1, quantile(avg_pred_0,.025), 
         1, quantile(avg_pred_0,.975), lwd = 2)
points(1, mean(avg_pred_0), pch = 16, cex = 1.2)
segments(2, quantile(avg_pred_1,.025), 
         2, quantile(avg_pred_1,.975), lwd = 2)
points(2, mean(avg_pred_1), pch = 16, cex = 1.2)
segments(3, quantile(avg_pred_2,.025), 
         3, quantile(avg_pred_2,.975), lwd = 2)
points(3, mean(avg_pred_2), pch = 16, cex = 1.2)
dev.off()

pdf("probabilities_delta_1.pdf", height = 4, width = 5) # Figure 6.1
par(mar=c(5,4,4,2))
plot("",
     xlim = c(-0.2,0.5),
     ylim = c(0,7),
     xlab = expression(paste(Delta," Average marginal probability")),
     ylab = "Density",
     main = "First difference: Scenario 2 - Scenario 1",
     axes = FALSE)
axis(1,at = seq(-.2,0.5, by = .1))
axis(2,at = seq(0,7, by = .5))
polygon(density(pred_delta_1), col = alpha("lightgrey", 0.85), border = NA)
segments(quantile(pred_delta_1,.025),0,quantile(pred_delta_1,.975),0, lwd = 2)
points(mean(pred_delta_1),0, pch = 16, cex = 1.2)
abline(v = 0, lty = 2, lwd = 1)
dev.off()

pdf("probabilities_delta_2.pdf", height = 4, width = 5) # Figure 6.2
par(mar=c(5,4,4,2))
plot("",
     xlim = c(-0.5,0.2),
     ylim = c(0,7),
     xlab = expression(paste(Delta," Average marginal probability")),
     ylab = "Density",
     main = "First difference: Scenario 3 - Scenario 2",
     axes = FALSE)
axis(1,at = seq(-.5,0.2, by = .1))
axis(2,at = seq(0,7, by = .5))
polygon(density(pred_delta_2), col = alpha("lightgrey", 0.85), border = NA)
segments(quantile(pred_delta_2,.025),0,quantile(pred_delta_2,.975),0, lwd = 2)
points(mean(pred_delta_2),0, pch = 16, cex = 1.2)
abline(v = 0, lty = 2, lwd = 1)
dev.off()


#================================================================
# D) Interaction effects

### Model 3
fit.3 <- stan_glmer(court_decision ~ auth_gov * obj_gov + auth_gov * obj_opp + cross_party + policy_area + (1 | law_id),
                    data = mydata, family = binomial(link = "logit"),
                    # priors (default November 2018)
                    # prior = normal(0,2.5), prior_intercept = normal(0,10),
                    chains=4, seed=123, iter=6000, warmup=1000)
print(fit.3)
#summary(fit.3)
# priors
prior_summary(fit.3)
# convergence
rhat.3 <- summary(fit.3)[, "Rhat"]
max(rhat.3)

# Sigma
var(unlist(ranef(fit.3)))
# Get full results for fixed effect coefficients
BETA.3 <- as.matrix(fit.3)[,2:length(fixef(fit.3))]
mean.BETA.3 <- apply(BETA.3,2,mean)
sd.BETA.3 <- apply(BETA.3,2,sd)
low.BETA.3 <- apply(BETA.3,2,quantile,c(0.025,0.975))[1,]
high.BETA.3 <- apply(BETA.3,2,quantile,c(0.025,0.975))[2,]
round(mean.BETA.3, digits = 2)
round(low.BETA.3, digits = 2)
round(high.BETA.3, digits = 2)


### Model 4
fit.4 <- stan_glmer(court_decision ~ z.days_passed * obj_gov + z.days_passed * obj_opp + cross_party + policy_area + (1 | law_id),
                    data = mydata, family = binomial(link = "logit"),
                    # priors (default November 2018)
                    # prior = normal(0,2.5), prior_intercept = normal(0,10),
                    chains=4, seed=123, iter=6000, warmup=1000)
print(fit.4)
#summary(fit.3)
# priors
prior_summary(fit.4)
# convergence
rhat.4 <- summary(fit.4)[, "Rhat"]
max(rhat.4)

# Sigma
var(unlist(ranef(fit.3)))
# Get full results for fixed effect coefficients
BETA.4 <- as.matrix(fit.4)[,2:length(fixef(fit.4))]
mean.BETA.4 <- apply(BETA.4,2,mean)
sd.BETA.4 <- apply(BETA.4,2,sd)
low.BETA.4 <- apply(BETA.4,2,quantile,c(0.025,0.975))[1,]
high.BETA.4 <- apply(BETA.4,2,quantile,c(0.025,0.975))[2,]
round(mean.BETA.4, digits = 2)
round(low.BETA.4, digits = 2)
round(high.BETA.4, digits = 2)

# Plot coefficients from Model 3
pdf("coefficients_3.pdf", height = 5, width = 10) # Figure 7
par(mar=c(5,16,3,3))
plot("", axes = FALSE,
     xlim = c(-5,5),
     ylim = c(0.5,10),
     xlab = "Posterior means with 95% HPDs",
     ylab = "", pch = 19, cex = 1.5,
     #main = "Coefficient posterior distributions with mean point estimates and 90% HPDs",,
     main = "GFCC strike")
axis(1,at = seq(-5,5,.5))
axis(2,at = seq(9,1,-2), las = 1, tick = TRUE,
     label = c("Author of policy in government",
               "Contested by governing lawmaker",
               "Contested by opposition lawmaker",
               "Author of policy in government x \n Contested by governing lawmaker",
               "Author of policy in government x \n Contested by opposition lawmaker"))
abline(h = seq(9,1,-2), lty = 2, lwd = .5, col = "grey")
# Add posterior densities
BETA.plot <- BETA.3[,c(1,2,3,14,15)]
for (i in 1:ncol(BETA.plot)) {
  BETA.density <- density(BETA.plot[,i])
  lower <- which(BETA.density$x < quantile(BETA.plot[,i], 0.00001))
  upper <- which(BETA.density$x > quantile(BETA.plot[,i], 0.99999))
  BETA.density$x <- BETA.density$x[-c(lower,upper)] # truncate density
  BETA.density$y <- BETA.density$y[-c(lower,upper)] # truncate density
  locate.density <- seq(9,1,-2)
  polygon(BETA.density$x,BETA.density$y+locate.density[i], col = alpha("lightgrey", 0.85), border=NA)
}
points(mean.BETA.3[c(1,2,3,14,15)],seq(9,1,-2), pch = 19, cex = 1.2)
segments(low.BETA.3[c(1,2,3,14,15)], seq(9,1,-2), high.BETA.3[c(1,2,3,14,15)], seq(9,1,-2), lwd = 2) # 95% HPDs
abline(v=0, lty = 2)
dev.off()


# Plot coefficients from Model 4
pdf("coefficients_4.pdf", height = 5, width = 10) # Figure 8
par(mar=c(5,16,3,3))
plot("", axes = FALSE,
     xlim = c(-5,5),
     ylim = c(0.5,10),
     xlab = "Posterior means with 95% HPDs",
     ylab = "", pch = 19, cex = 1.5,
     #main = "Coefficient posterior distributions with mean point estimates and 90% HPDs",,
     main = "GFCC strike")
axis(1,at = seq(-5,5,.5))
axis(2,at = seq(9,1,-2), las = 1, tick = TRUE,
     label = c("Days passed since adoption of policy",
               "Contested by governing lawmaker",
               "Contested by opposition lawmaker",
               "Days passed since adoption of policy x \n Contested by governing lawmaker",
               "Days passed since adoption of policy x \n Contested by opposition lawmaker"))
abline(h = seq(9,1,-2), lty = 2, lwd = .5, col = "grey")
# Add posterior densities
BETA.plot <- BETA.4[,c(1,2,3,14,15)]
for (i in 1:ncol(BETA.plot)) {
  BETA.density <- density(BETA.plot[,i])
  lower <- which(BETA.density$x < quantile(BETA.plot[,i], 0.00001))
  upper <- which(BETA.density$x > quantile(BETA.plot[,i], 0.99999))
  BETA.density$x <- BETA.density$x[-c(lower,upper)] # truncate density
  BETA.density$y <- BETA.density$y[-c(lower,upper)] # truncate density
  locate.density <- seq(9,1,-2)
  polygon(BETA.density$x,BETA.density$y+locate.density[i], col = alpha("lightgrey", 0.85), border=NA)
}
points(mean.BETA.4[c(1,2,3,14,15)],seq(9,1,-2), pch = 19, cex = 1.2)
segments(low.BETA.4[c(1,2,3,14,15)], seq(9,1,-2), high.BETA.4[c(1,2,3,14,15)], seq(9,1,-2), lwd = 2) # 95% HPDs
abline(v=0, lty = 2)
dev.off()


#================================================================
# E) Robustness checks

### Model 5
fit.5 <- stan_glmer(court_decision ~ senior_gov_cat + obj_opp + z.days_passed + cross_party + policy_area + (1 | law_id),
                    data = mydata, family = binomial(link = "logit"),
                    # priors (default November 2018)
                    # prior = normal(0,2.5), prior_intercept = normal(0,10),
                    chains=4, seed=123, iter=6000, warmup=1000)
print(fit.5)
summary(fit.5)
# priors
prior_summary(fit.5)
# convergence
rhat.5 <- summary(fit.5)[, "Rhat"]
max(rhat.5)

# Sigma
var(unlist(ranef(fit.5)))
# Get full results for fixed effect coefficients
BETA.5 <- as.matrix(fit.5)[,2:length(fixef(fit.5))]
mean.BETA.5 <- apply(BETA.5,2,mean)
sd.BETA.5 <- apply(BETA.5,2,sd)
low.BETA.5 <- apply(BETA.5,2,quantile,c(0.025,0.975))[1,]
high.BETA.5 <- apply(BETA.5,2,quantile,c(0.025,0.975))[2,]
round(mean.BETA.5, digits = 2)
round(low.BETA.5, digits = 2)
round(high.BETA.5, digits = 2)


### Model 6
fit.6 <- stan_glmer(court_decision ~ z.obj_gov_time_as_mp + obj_opp + z.days_passed + cross_party + policy_area + (1 | law_id),
                    data = mydata, family = binomial(link = "logit"),
                    # priors (default November 2018)
                    # prior = normal(0,2.5), prior_intercept = normal(0,10),
                    chains=4, seed=123, iter=6000, warmup=1000)
print(fit.6)
summary(fit.6)
# priors
prior_summary(fit.6)
# convergence
rhat.6 <- summary(fit.6)[, "Rhat"]
max(rhat.6)

# Sigma
var(unlist(ranef(fit.6)))
# Get full results for fixed effect coefficients
BETA.6 <- as.matrix(fit.6)[,2:length(fixef(fit.6))]
mean.BETA.6 <- apply(BETA.6,2,mean)
sd.BETA.6 <- apply(BETA.6,2,sd)
low.BETA.6 <- apply(BETA.6,2,quantile,c(0.025,0.975))[1,]
high.BETA.6 <- apply(BETA.6,2,quantile,c(0.025,0.975))[2,]
round(mean.BETA.6, digits = 2)
round(low.BETA.6, digits = 2)
round(high.BETA.6, digits = 2)


# Predicted probabilities
fixef(fit.6)
tmpdata <- as.data.frame(cbind(1, mydata$z.obj_gov_time_as_mp, mydata$obj_opp, mydata$z.days_passed, mydata$cross_party,
                               ifelse(mydata$policy_area == "education/research", 1,0),
                               ifelse(mydata$policy_area == "family", 1,0), ifelse(mydata$policy_area == "healthcare", 1,0), ifelse(mydata$policy_area == "interior", 1,0),
                               ifelse(mydata$policy_area == "labour/social insurance", 1,0), ifelse(mydata$policy_area == "other", 1,0), ifelse(mydata$policy_area == "public finances", 1,0),
                               ifelse(mydata$policy_area == "rights", 1,0), ifelse(mydata$policy_area == "transport/public infrastructure", 1,0)))
dim(tmpdata)
colnames(tmpdata) <- names(fixef(fit.6))

# Set up design matrix Z for random effects
N <- dim(mydata)[1]
q <- dim(ranef(fit.6)[[1]])[1]
Z <- as.data.frame(matrix(NA, nrow = N, ncol = q))
dim(Z)
colnames(Z) <- rownames(ranef(fit.6)[[1]])
for (i in colnames(Z)) {
  Z[,i] <- ifelse(mydata$law_id == i,1,0) 
}

# Check whether data matrix and design matrix are set up correctly
mydata[1,c("case","law_id", "z.obj_gov_time_as_mp","obj_opp","cross_party","policy_area")]
tmpdata[1,] # values for first observation should match mydata
Z[1,as.character(mydata[1,"law_id"])] # should return 1 if correct

# Fixed effect coefficients from all sampling iterations
BETA <- as.matrix(fit.6)[,1:length(fixef(fit.6))]
dim(BETA)
BETA[1,]

# Random effects from all sampling iterations
start.alpha <- length(fixef(fit.6)) + 1
GAMMA <- as.matrix(fit.6)[,start.alpha:dim(as.matrix(fit.6))[2]]
GAMMA <- GAMMA[,-dim(GAMMA)[2]] # drop sigma
dim(GAMMA)


# I should set up a function for the following
# Set up data matrices
summary(mydata$z.obj_gov_time_as_mp)
X.0 <- as.matrix(tmpdata)
X.0[,match("z.obj_gov_time_as_mp",names(fixef(fit.6)))] <- -0.12

X.1 <- as.matrix(tmpdata)
X.1[,match("z.obj_gov_time_as_mp",names(fixef(fit.6)))] <- 0.5

X.2 <- as.matrix(tmpdata)
X.2[,match("z.obj_gov_time_as_mp",names(fixef(fit.6)))] <- 1.0

X.3 <- as.matrix(tmpdata)
X.3[,match("z.obj_gov_time_as_mp",names(fixef(fit.6)))] <- 1.5

X.4 <- as.matrix(tmpdata)
X.4[,match("z.obj_gov_time_as_mp",names(fixef(fit.6)))] <- 2.0

X.5 <- as.matrix(tmpdata)
X.5[,match("z.obj_gov_time_as_mp",names(fixef(fit.6)))] <- 2.5

X.6 <- as.matrix(tmpdata)
X.6[,match("z.obj_gov_time_as_mp",names(fixef(fit.6)))] <- 3.0

Z <- as.matrix(Z)

linpred_0 <- X.0 %*% t(BETA) + Z %*% t(GAMMA) # Linear predictor with X.0
dim(linpred_0)
pred_0 <- exp(linpred_0)/(1+exp(linpred_0)) # Through link function
dim(pred_0)
avg_pred_0 <- apply(pred_0,2,mean)
mean(avg_pred_0)
quantile(avg_pred_0,c(0.025,0.975))

linpred_1 <- X.1 %*% t(BETA) + Z %*% t(GAMMA) # Linear predictor with X.1
dim(linpred_1)
pred_1 <- exp(linpred_1)/(1+exp(linpred_1)) # Through link function
dim(pred_1)
avg_pred_1 <- apply(pred_1,2,mean)
mean(avg_pred_1)
quantile(avg_pred_1,c(0.025,0.975))

linpred_2 <- X.2 %*% t(BETA) + Z %*% t(GAMMA) # Linear predictor with X.2
dim(linpred_2)
pred_2 <- exp(linpred_2)/(1+exp(linpred_2)) # Through link function
dim(pred_2)
avg_pred_2 <- apply(pred_2,2,mean)
mean(avg_pred_2)
quantile(avg_pred_2,c(0.025,0.975))

linpred_3 <- X.3 %*% t(BETA) + Z %*% t(GAMMA) # Linear predictor with X.3
dim(linpred_3)
pred_3 <- exp(linpred_3)/(1+exp(linpred_3)) # Through link function
dim(pred_3)
avg_pred_3 <- apply(pred_3,2,mean)
mean(avg_pred_3)
quantile(avg_pred_3,c(0.025,0.975))

linpred_4 <- X.4 %*% t(BETA) + Z %*% t(GAMMA) # Linear predictor with X.3
dim(linpred_4)
pred_4 <- exp(linpred_4)/(1+exp(linpred_4)) # Through link function
dim(pred_4)
avg_pred_4 <- apply(pred_4,2,mean)
mean(avg_pred_4)
quantile(avg_pred_4,c(0.025,0.975))

linpred_5 <- X.5 %*% t(BETA) + Z %*% t(GAMMA) # Linear predictor with X.3
dim(linpred_5)
pred_5 <- exp(linpred_5)/(1+exp(linpred_5)) # Through link function
dim(pred_5)
avg_pred_5 <- apply(pred_5,2,mean)
mean(avg_pred_5)
quantile(avg_pred_5,c(0.025,0.975))

linpred_6 <- X.6 %*% t(BETA) + Z %*% t(GAMMA) # Linear predictor with X.3
dim(linpred_6)
pred_6 <- exp(linpred_6)/(1+exp(linpred_6)) # Through link function
dim(pred_6)
avg_pred_6 <- apply(pred_6,2,mean)
mean(avg_pred_6)
quantile(avg_pred_6,c(0.025,0.975))


# Plot predicted probabilities
pdf("jitter_plot.pdf", height = 5, width = 10) # Figure 5 Appendix
set.seed(123)
plot(jitter(mydata$z.obj_gov_time_as_mp, 2), jitter(mydata$court_decision, 0.05), pch = 16,
     xlab = "Time as Member of Parliament (standardized)",
     ylab = "Predicted probability",
     main = "Predicted probabilities of Strike = 1")
points(c(-0.12,0.5,1.0,1.5,2.0,2.5,3.0),
       c(mean(avg_pred_0),mean(avg_pred_1),mean(avg_pred_2),
         mean(avg_pred_3),mean(avg_pred_4),mean(avg_pred_5),mean(avg_pred_6)), pch = 19, cex = 1.2)
segments(c(-0.12,0.5,1.0,1.5,2.0,2.5,3.0),
         c(quantile(avg_pred_0,0.025),quantile(avg_pred_1,0.025),quantile(avg_pred_2,0.025),
           quantile(avg_pred_3,0.025),quantile(avg_pred_4,0.025),quantile(avg_pred_5,0.025),quantile(avg_pred_6,0.025)), 
         c(-0.12,0.5,1.0,1.5,2.0,2.5,3.0),
         c(quantile(avg_pred_0,0.975),quantile(avg_pred_1,0.975),quantile(avg_pred_2,0.975),
           quantile(avg_pred_3,0.975),quantile(avg_pred_4,0.975),quantile(avg_pred_5,0.975),quantile(avg_pred_6,0.975)), lwd = 2) # 95% HPDs
dev.off()



### Model 7
table(mydata$proceed)
mydata$proceed[mydata$proceed == "organ"] <- "abstract"
mydata$proceed <- as.factor(mydata$proceed)
mydata$proceed <- relevel(mydata$proceed, ref = "const co")

fit.7 <- stan_glmer(court_decision ~ obj_gov + obj_opp + proceed + (1 | law_id),
                    data = mydata, family = binomial(link = "logit"),
                    # priors (default November 2018)
                    # prior = normal(0,2.5), prior_intercept = normal(0,10),
                    chains=4, seed=123, iter=6000, warmup=1000)
print(fit.7)
summary(fit.7)
# priors
prior_summary(fit.7)
# convergence
rhat.7 <- summary(fit.7)[, "Rhat"]
max(rhat.7)

# Sigma
var(unlist(ranef(fit.7)))
# Get full results for fixed effect coefficients
BETA.7 <- as.matrix(fit.7)[,2:length(fixef(fit.7))]
mean.BETA.7 <- apply(BETA.7,2,mean)
sd.BETA.7 <- apply(BETA.7,2,sd)
low.BETA.7 <- apply(BETA.7,2,quantile,c(0.025,0.975))[1,]
high.BETA.7 <- apply(BETA.7,2,quantile,c(0.025,0.975))[2,]

pdf("robustness.pdf", height = 6, width = 10) # Figure 6 Appendix
par(mar=c(5,14,3,3))
plot("", axes = FALSE,
     xlim = c(-3,2),
     ylim = c(0.5,9.5),
     xlab = "Posterior means with 95% HPDs",
     ylab = "", pch = 19, cex = 1.5,
     #main = "Coefficient posterior distributions with mean point estimates and 90% HPDs",,
     main = "GFCC strike")
axis(1,at = seq(-3,2,.5))
axis(2,at = seq(7,1,-2), las = 1, tick = TRUE,
     label = c("Contested by governing lawmaker",
               "Contested by opposition lawmaker",
               "Abstract review",
               "Concrete review"))
abline(h = seq(7,1,-2), lty = 2, lwd = .5, col = "grey")
# Add posterior densities
for (i in 1:length(colnames(BETA.7))) {
  BETA.density <- density(BETA.7[,i])
  lower <- which(BETA.density$x < quantile(BETA.7[,i], 0.00001))
  upper <- which(BETA.density$x > quantile(BETA.7[,i], 0.99999))
  BETA.density$x <- BETA.density$x[-c(lower,upper)] # truncate density
  BETA.density$y <- BETA.density$y[-c(lower,upper)] # truncate density
  locate.density <- seq(7,1,-2)
  polygon(BETA.density$x,BETA.density$y+locate.density[i], col = alpha("lightgrey", 0.85), border=NA)
}
points(mean.BETA.7,seq(7,1,-2), pch = 19, cex = 1.2)
segments(low.BETA.7, seq(7,1,-2), high.BETA.7, seq(7,1,-2), lwd = 2) # 95% HPDs
abline(v=0, lty = 2)
dev.off()

